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In this article we study the ferromagnetic behavior of ABC-stacked trilayer graphene. This is done 
using a nearest-neighbor tight-binding model, in the presence of long-range Coulomb interactions. 
For a given electron-electron interaction g and doping level n, we determine whether the total 
energy is minimized for a paramagnetic or ferromagnetic configuration of our variational parameters. 
The g versus n phase diagram is first calculated for the unscreened case. We then include the 
effects of screening using a simplified expression for the fermion bubble diagram. We show that 
ferromagnetism in ABC-trilayer graphene is more robust than in monolayer, in bilayer, and in 
ABA-trilayer graphene. Although the screening reduces the ferromagnetic regime in ABC-trilayer 
graphene, the critical doping level remains one order of magnitude larger than in unscreened bilayer 
graphene. 

PACS numbers: 75.70.Cn, 73.22.Pr, 73.20.At 



I. INTRODUCTION 

Within a decade after the discovery of graphene flakes 
by mechanical exfoliation,^ numerous methods have been 
developed to create larger and cleaner samples, realized 
both as single layers and as stacked layers of grapheneP^ 

Early on, it was realized that stacked graphene lay- 
ers behave differently than both a single layer and 3D 
graphite. For example, in bilayer graphene the disper- 
sion is quadratic instead of linear and the electrons be- 
have as massive chiral particles, which is a completely 
new type of particle. Few-layer graphene is still a 2D 
system, hence the quantum Hall effect can be observed. 
For monolayer graphene, the plateaus in the Hall con- 
ductivity are located at half integer multiples of 4e 2 //i, 8 
originating from a Landau level at zero energy which is 
half filled by electrons and half filled by holes. In bilayer 
graphene, this particular Landau level has an extra de- 
generacy resulting in Hall plateaus at integer values of 
4e 2 /h and a quantum Hall effect that is different from 
the one in a monolayer as well as from the quantum Hall 
effect found in usual two dimensional electron gases P In 
addition to the number of layers, the order of the stacking 
also influences the physical properties significantly. 

In multilayer graphene, the different layers can have 
three distinct orientations with respect to the bottom 
one. Bernal stacking (or AB stacking) is the configura- 
tion in which the B sublattice of the odd layers are oppo- 
site to the A sublattice of the even layers. The Hamilto- 
nian of a system with an even number 2N of layers can 
be rewritten in a block diagonal form, where the N differ- 
ent blocks are bilayer-like Hamiltonians. The blocks can 
be linked by hopping parameters that couple lattice sites 
on next-nearest planes. For an odd number (27V + 1) of 
layers, one of the blocks is the monolayer Hamiltonian. 
Therefore, these systems have a linear band in addition 
to the N parabolic onesP^ 

In ABC stacked multilayer graphene, the B sublattice 
of each layer lies opposite to the A sublattice of the layer 




Figure 1. (Color online) Atomic structure of Ai?C-trilayer 
graphene. 



above it, but opposite to the honeycomb centers in the 
layer beneath it (see Fig. [I]). Since electrons that are 
placed oppositely in two bordering planes dimerize, re- 
sulting in an energy shift away from zero, these multilay- 
ers can, for low energies, be described by a 2 x 2 effective 
matrix Hamiltonian, which is governed by the indirect 
(effective) hopping between the two atoms in the outer 
planes that have no neighbor in the adjacent layer. This 
effective hopping is a process consisting of N — 1 inter- 
plane nearest-neighbor hoppings, combined with N in- 
plane nearest-neighbor hoppings, resulting in an energy 
dispersion around the K-points, En ~ v^k N /t^~^^ 
A tight-binding approach for an increasing number of 
layers should in principle include hopping between more 
distant carbon atoms. The longknown Slonczewski- 
Weiss-McClure (SWMc) modeP^U^ accounts for next- 
nearest-neighbor hopping, as well as hopping between 
next-nearest planes. In fact, trilayer graphene can be 
used to obtain the values of the different hopping param- 
eters by fitting experimental data to the SWMc modelJ^ 
However, often it is sufficient to take into account only 
the intra- and interplane nearest-neighbor hopping pa- 
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rameters. 

Recent experimental and theoretical studies of trilayer 
graphene have shown that magnetotransport and elec- 
tronic transport properties,^ thermoelectric transport 
propertied, and chiral tunneling^ indeed depend on 
the stacking order. Furthermore, one can open a sizeable 
bandgap in AB C-st&cked trilayers (120 meV) by apply- 
ing an external electric field, while for an AB A-tiiidiyei 
no gap is observed under the same conditions.^ 

Extensive research into the band structure of ABC- 
multilayer graphene has been done recently using an ef- 
fective mass approximation.^ It was found that the elec- 
tron and hole bands touching at zero energy support 
chiral quasiparticles characterized by a Berry phase of 
Ntt for N layers. The phonon spectrum of A£?C-stacked 
graphene has been investigated theoretically using den- 
sity functional theory^ anc [ experimentally by using in- 
frared absorption spectroscopy, where the intensities have 
been found to be much stronger than that of bilayer 
graphene. 2 2 Using magnetic fields up to 60T, there has 
been evidence of the integer quantum Hall effect in tri- 
layer graphene. 22 The Hall resistivity plateaus have been 
reproduced by using a self-consistent Hartree calculation 
on A£>C-stacked graphene.^! It has been suggested that 
the differences in the quantum Hall effect between ABC- 
and ABA- stacking might be used to identify the stacking 
order of high-quality trilayer samples™ By using infrared 
absorption spectroscopy, it has been shown that the op- 
tical conductivity spectra for ABC- and ABA- stacked 
graphene differs considerably.^ These optical properties 
have been calculated and reproduced in the framework 
of a tight-binding modelP^ Finally, it can be mentioned 
that high-resolution transmission microscopy of ABC- 
stacked trilayer graphene on a SiC surface has success- 
fully provided information on the interlayer distances of 
ABC- trilayer grapheneP^ 

In this article we investigate the magnetic properties of 
ABC- trilayer graphene by using a nearest-neighbor tight- 
binding model, in the presence of long-range Coulomb 
interactions. For interacting electrons, the system can 
gain energy by aligning the spins of the electrons. This 
exchange mechanism is accompanied by a cost in kinetic 
energy due to the Pauli principle. After fixing the dop- 
ing level and interaction strength, one can use a varia- 
tional approach to determine whether the system sponta- 
neously magnetizes or remains paramagnetic. For mono- 
layer graphene, the system only magnetizes if the inter- 
action strength is tuned to unphysically high values. De- 
pending on the doping level n, this phase transition can 
be first or second order™ For bilayer graphene the sys- 
tem can be ferromagnetic for the estimated value of the 
Coulomb interaction (g = 2.1), but the electron den- 
sity has to be as low as n ~ 10 9 cm -2 for the material 
to become magnetic.^ This is on the brink of what is 
experimentally achievable, since it is not possible to cre- 
ate perfectly undoped graphene in experiment, due to 
the formation of electron hole puddles^ and impurities 
trapped in the substrate. In AB A- trilayer, the inter- 
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Figure 2. (Color online) (a) Numerically calculated, full dis- 
persion of Ai^C-trilayer graphene (two lowest bands), (b) 
Zoom-in on one of the K-points. (c) Numerically calculated 
low energy approximation of ABC-trii&yer graphene disper- 
sion (expansion around the K-point). 



play between the linear and the parabolic bands opens 
up possibilities for both spin,- and band-ferromagnetism, 
but only at low electron doping.^ 

Although in a low energy approximation A£?C-trilayer 
graphene seems to be - in a way - the three layer gener- 
alization of the Bernal stacked bilayer™ it is worth a fur- 
ther investigation because its cubic energy dispersion is 
expected to enhance significantly the phase-space where 
the ferromagnetic regime occurs. In addition, screening 
should play an important role, due to the diverging den- 
sity of states. Here we show that this is indeed the case: 
although the screening reduces the regime of parameters 
for the occurrence of ferromagnetism, the latter remains 
at least one order of magnitude more robust than in un- 
screened bilayer graphene. The outline of our paper is 
the following: we set up the model in Sec.[Hj present our 
results of the unscreened case in Sec. IIIIl and look at the 



effects of screening in Sec. IV Our conclusions are drawn 
in Sec. N\ 



II. THE MODEL 

We use a tight-binding model which takes into account 
the hopping of electrons to nearest-neighbor inplane and 
interplane sites. In real space, the Hamiltonian is given 

by 



H = Hq + Hj , 
with the non-interacting part being 
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Figure 3. (Color online) Sketch of the electron- (hole-) pockets for three configurations of the system - (a) paramagnetic, (b) 
ferromagnetic with one type of charge carrier and (c) ferromagnetic with two types of carriers. 



where i and j label the lattice sites, a G {t, 1} labels spin, 
n G {1, 2, 3} labels the layer, t « 3 eV denotes the intra- 
layer nearest-neighbor hopping parameter, t± « 0.35 eV 
denotes the interlayer nearest-neighbor hopping, and the 
operator (c) creates (annihilates) an electron on sub- 
lattice C G {A, B}. Hj is the interaction Hamiltonian. 
Since the stacking considered is ABC the A sublattice 
in the bottom layer (layer 1) and the B sublattice in the 
top layer (layer 3) do not have direct neighbors in an 
adjacent layer. The electrons interact via a Coulomb in- 
teraction, which can be included in our model by the 
term 

Hi = \j d 2 x C Z 2 y{^ D (x-y)bi(x)p 1 (y) + p 2 (x) /02 (y) 

+ P 3 (x)p 3 (y)] + ^ ND (x - y)[pi(x) P2 (y) 

+ P2(x)pi(y) + p 2 (x)p 3 (y) + Ps(x)p 2 (y)] 

+ ^ 2ND (x - y)[p 1 (x) P3 (y) + P 3 (x)pi(y)]}, (3) 

where the density of electrons in the n-th layer is 
given by p n (x) = ^^ n (x)i )n (x), with ^, n (x) = 
(a C7)n (x),6 . >n (x)), where a a , n (x) and 6 a , n (x) are the 
field operators corresponding to di i(J ^ n and 6i ><T>n , respec- 
tively. The interaction potentials for the in-plane (D), 
the nearest-neighbor planes (ND) and the next-nearest- 
neighbor planes (2ND) are given by 



tonian acquires the form 



V D (x-y) 
^ ND (x-y) 
y2ND (x _ y) 
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In these interaction potentials, d ~ 3.2 A is the inter- 
layer distance, e the electron charge, and e the dielectric 
constant of the substrate. 



A. Kinetic energy 

After Fourier transforming and expanding the mo- 
menta around the If -point, the non-interacting Hamil- 



dk*t(k)W(k)* a (k), (4) 

where G creates a particle with momentum k on sub- 
lattice c G {a, b} in layer n with spin a, and H is a 6 x 6 
matrix given by 
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where u = ke % ^\ In the above expression, k = |k| 
is the norm of the two-dimensional momentum vector, 
0(k) = arctan (k y /k x ) is the angle of the momentum vec- 
tor, hvp = (3/ 2) at is the Fermi velocity in terms of the 
lattice constant a = 1.42 A and intralayer hopping pa- 
rameter t, and 7i = t±/(hvF)- 

Although it is possible to write an analytic expression 
for the low energy approximation of the single-particle 
dispersion for A£?C-trilayer graphenep^ this is not the 
case for the required diagonalization matrix for T~L. For 
this reason, we calculate both numerically. The full dis- 
persion is shown in Fig. [2|a) and (b), together with an 
expansion of the energy bands around the K-point (i.e. 
eigenvalues of Eq. (|5|), which are indeed cubic for small 
momenta (at small momenta E(k) « ^z(vp/t']_)k 3 for the 
two lowest bands), see Fig. [2jc). 

When the system undergoes a phase transition into a 
ferromagnetic state, pockets of one spin configuration - 
let us say up - will be larger than the pocket of spin-down 
electrons [see Fig. |3ja)-(b)] Moreover, it is also possible 
to have two types of charge carriers in the system, i.e. 
the formation of spin-up electron-pockets and spin-down 
hole pockets [see Fig. |3jc)]. 

To compute the energy of an electron or hole pocket 
of size Q a (see Fig. |3|, we have to compute the integral 

rE{Q a ) 

AK = / EV(E)dE, 
Jo 
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where V(E) is the density of states 

with A denoting the area of the unit cell and N is the 
number of states below E. We compute the inverse of the 
dispersion relation E{k) numerically. Note that for small 
pocket sizes, AK ~ Q^. When compared with mono- 
layer graphene (AK m,L ~ Q&r^ and bilayer graphene 
(AK b - L ~ Q^jM[t is evident that the kinetic energy cost 
of an electron (hole) pocket is smaller in A£?C-trilayer 
graphene than in the fewer-layered carbon structures. 




-x' [10- 4 ] 



B. Exchange energy 

When calculating the energy contribution coming from 
Hi, the direct contribution (i.e. the Hartree term) can- 
cels due to the positive Jellium background. The only 
term left is the exchange contribution (i.e. the Fock 
term), which favors spin alignment. However, spin align- 
ment will result in a cost in kinetic energy due to the 
Pauli exclusion principle. Thus, ferromagnetism will oc- 
cur or not, depending on the competition between the 
kinetic energy and the exchange energy. 

In the Appendix, it is shown that the exchange energy 
of a configuration as in Fig. [3j where the spin- up and the 
spin-down bands fill up differently, can be written in a 
way similar to the one in bilayer graphene, 28 



A 



dk dk f 

(2^)2 (2^p 



EE E 

cr,a s=l a, (3=1 
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x^(k',k) x i a (k,k , )K(k' - kH, a)0 (k'K i/3)a (k) 



Here, a and /3 label the band index and a labels the valley, 
but we will neglect intervalley scattering and only focus 
on the K point. 7V,a,a(k) are the Fermi functions and 
the expressions for V s (k 7 — k) are given in the Appendix. 
In comparison with the bilayer, there are six x matrices 
instead of two and they are no longer 4x4, but 6x6. 
Moreover, they can only be computed numerically (see 
the Appendix for more details). 

Since we have exp anded around the K point, we intro- 
duce a cutoff A = ^2tt/A in such a way that the num- 
ber of states in the Brillouin zone is conserved. Using 
the cutoff, we can measure momenta (and hence pocket 
sizes) in units of A and energies in units of HvfA(~ 7.2 
eV) . This makes all our variables and parameters dimen- 
sionless and after setting h = 1, vp = 1, and A = 1 they 
have the following values: t = 0.42, t± = 0.05, a = 1.56, 
andd = 3.7P 



Figure 4. (Color online) AE versus — x at electron-electron 
coupling g = 6. The phase transition occurs at the dop- 
ing Qd ~ 0.0116 that produces the thick black curve. Inset: 
Minimum energy versus doping Qd. The phase transition is 
identified by the value of Qd where AE m in first becomes non- 
zero. 



III. UNSCREENED CASE 

A. Numerical solution 

The exchange energy E ex /A given by Eq. (|6| is 
solved numerically using the double exponential (DE) 
algorithm^ (the DE algorithm is originally intended for 
ID integrals, but is extended to 3D to perform the ex- 
change integrals). Due to the singular behavior of the 
Coulomb potentials, the integral must undergo a series 
of transformations. Firstly, the integral is transformed 
to polar coordinates, where we introduce a cutoff A for 
integrals over the norm of the the momentum. A change 
of variables is then applied, such that these integrations 
range from zero to one. This permits the singular be- 
havior along k = k 7 to be rotated by a Duffy coordinate 
transformational 



F(k,k',6) 



2tt pi pi 

d0 dk dk'- , 

o Jo Jo \J (k /2 — Ikk'Q cos 9 + Q 2 k 2 

2tt r l pi 



= / d6 dk dk" 
Jo Jo Jo 

F(kk",k,9) 



F(k,kk",9) 



y/l - 2k"Qcos9 + Q 2 k 



\/k"' 2 -2k"Q cos + Q 2 
(7) 



This formula is derived by splitting the k' integration 
into two separate integrations from to k and from k 
to 1. Making the change of variables k' = kk" on the 
first integral leads to the first term on the right hand 
side of Eq. In the second integral, with integration 
boundaries k and 1, the identity Jq dk dk' f(k,k') = 

Jq 1 dk Jq dk' f(k' \k) is applied. Thus, a change of vari- 
ables k' — kk" leads to the second term on the right 
hand side of Eq. ([F]). 
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The singularities are now confined to lines parallel to 
the /c-axis. However, there are now two such lines of 
singularities in the integrand, located at k" = h\ ^ 1 
and k" = Y12 7^ 1. The lines of singularities located at h\ 
and /12 must be moved to k" = 1 by a change of variables. 
After the change of variables, the integration boundaries 
are no longer confined to zero and one. Since the DE 
algorithm is only capable of handling singularities at the 
integration boundaries, all integrals are split at k" = 1 
(where the singularities are now located), before being 
performed. 

The Hamiltonian matrix T~L of Eq. ([5| is diagonalized 
numerically using the Jacobi diagonalization algorithm, 
which is extended to handle a Hermitian 6x6 matrix 
by solving the corresponding 12 x 12 real symmetric 
matrix.^ The resulting diagonalization matrix M. (k) is 
used inside the \ matrices of Eq. Q to calculate the ex- 
change energy, while the resulting dispersion E(k) is used 
to calculate the kinetic energy (see Appendix for details) . 

The numerical diagonalization process does not pro- 
vide E~ 1 (k), which is needed to calculate the kinetic en- 
ergy. Thus, the inverse is approximated by linear inter- 
polation of the dispersion. Integration by parts yields 
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AK = E(Q)N(E(Q)) 



E(Q) 



N(E)dE, 



which is used in order to avoid explicit numerical evalu- 
ation of dN/dE. 

Consider a paramagnetic state with doping Qd and a 
ferromagnetic state with electron (or hole) pockets 
and Then, the kinetic energy difference is calculated 

by 



AE k 



A 



[AK(Qf) + AK(Qi) - 2AK(Q d )] . 



The difference in exchange energy AE ex /A is calculated 
by subtracting E ex /A of the paramagnetic state from the 
corresponding energy of the ferromagnetic state. For an 
unperturbed system, both spin channels are filled up to 
the Fermi- momentum Qd [see Fig. j3^a)]. Due to the ex- 
change mechanism, the system can prefer a ferromagnetic 
state with either one type of carrier or two types of car- 
riers [see Fig. [3jb)-(c)]. These perturbations are param- 
eterized by the variable x, which is positive for one type 
of carrier and given by 

Q^ = d Q^ x. 

For two types of carriers, x is defined to be negative and 
parameterizes the electron and hole pocket as 

Q*=2Q 2 d + \x\, Ql = \x\, 

where we assume the electron pocket in the spin-up chan- 
nel. Using this parametrization, particle conservation 
is satisfied. It is convenient to introduce x' = x — Q 2 d1 
such that x' = represents the unperturbed state (i.e. 
Qt = Qi = Qd)- Then, AE/A = AE kin /A + AE ex /A 
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Figure 5. (Color online) Phase diagram for ABC-trilayer 
graphene, in the case of an unscreened potential (solid blue 
line) and a screened potential (dashed black line). The red 
dots are the calculated values, while the solid line is an inter- 
polation function based on the calculated points. 



can be plotted as a function of x' for given electron- 
electron coupling g and doping Qd (see Fig. 
The minimum of AE{x)/A is estimated numerically by 
interpolation of points close to the minimum. The criti- 
cal doping, where the minimum AE m i n /A of AE(x)/A is 
zero, is found numerically by solving AE min (Qd)/A = 0. 
Since each minimum is a time consuming calculation, a 
simple binary search pattern is used (see inset of Fig. [4| . 



B. Phase diagram 

For a fixed value of g = 6, we see in Fig. [4] the behav- 
ior of AE as a function of pocket sizes, upon varying the 
doping Qd. For some doping values, AE is positive defi- 
nite (paramagnetic phase), while for others AE attains a 
negative minimum (ferromagnetic phase). Inspection of 
the critical curve (thick line) shows that there is a first 
order phase transition between the paramagnetic and fer- 
romagnetic phases. Repeating the entire procedure for 
different values of g leads to the g versus n phase diagram 
depicted in Fig.pl where n = Qj/2. The continuous solid 
line is an interpolation function of the calculated points. 

These results were obtained by neglecting higher order 
corrections that lead to screening of the Coulomb poten- 
tial. These effects will be considered in the next section. 



IV. EFFECTS OF SCREENING 

A. Screened potential 

Fourier transforming the real-space potentials V D , 
V ND and V 2ND and going to dimensionless variables 
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B. Phase diagram 

Notice that Eq. ([9| does not converge to the true un- 
screened potential V mn as U r i — ^ 0. In order to achieve 
such a convergence, Eq. ([9| must be changed to 



V mn (k) 



V{k)li tot {\L,iuY 



where 



TL tot {k,ko) = ^n(Q ff) k,iu;) = 5^II r ,(Q (T ,k ) iw). 



Figure 6. (Color online) Plot of the bilayer graphene polariza- 
tion II(A;f, k, 0). The analytical expression derived in Ref. 1371 
is shown as a blue dashed line for Jzf — 0.025 and as a red 
dotted line for — 0.05. The black solid line is the linear 
asymptote nk valid at large k, which we here extrapolate to 
small k. 



yields 
V D — 



d _ 27r 9 



rND _ 2 ^ge kd 



k ' 



ylMJ = 



k 



V 



2ND 



2irge 



-2kd 
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where g = e 2 /ehvF- As can be seen from Eq. (A2) in 
the Appendix, the bare interaction line of A£?C-trilayer 
graphene becomes a matrix F mn , where m and n are 
layer indices. Therefore, the RPA renormalization of the 
potential can be described by the Dyson-like equation 



m n 



(8) 



where r and I are layer indices. Let V mn be the renor- 
malized potential. Then, 

^^y\^Q^\ = - V mr (k)U r i(k, z<j)Vj n (k, iu), 

rl m T 1 rl 

where 

xGo'V-k,^-^) 

and Gq ,1t is the non-interacting Green's function of the 
system. Eq. ([8| is difficult to solve due to the layer depen- 
dence. However, for sufficiently low momenta e~ kd ~ 1 
and e~ 2kd ~ 1, which means that Vij ~ V = 2ng/k. 
Thus, the layer dependence is removed, and Eq. (|8| can 
be solved with respect to Vij = V: 



V(k) 



V(k) 



i-v(k) E P «n t . I (k ) iw)' 



(9) 



Since we are only interested in the long wavelength be- 
havior, then uj — )> 0. For both, monolayer and bilayer 
graphene, the polarization n(Q cr ,k, 0) behaves linearly 
in k for large /c, independent of Fermi momentum Q a , 
and exhibits an identical sloped This occurs because 
the dispersions are linear in the large-/c limit for both 
systems, and the Green's functions depend on the dis- 
persion. Since the dispersion of A£?C-trilayer graphene 
is also linear in the large- k limit with the same slope 
as of the single- and bilayer dispersions, it is reason- 
able to assume that the linear behavior of n(Q cr , k, 0) is 
also present for A£?C-trilayer graphene. In the exchange 
energy integrations, there are terms that are integrated 
from zero to the edge of the Brillouin zone (i.e. the cutoff 
A = 1). Therefore, we will first focus on the screening ef- 
fects coming from the linear behavior of U(Q a , k, 0) and 
approximate it by IT^Q^k, 0) = nk. An analytical ex- 
pression of n(Q cr , k, 0) was calculated by GamayurP^for 
bilayer graphene and is plotted in Fig. [6] for two values 
of the Fermi momentum (dashed and dotted lines) and 
compared with the linear estimate, where n ~ —0.12495 
(solid line) Notice that the high-A: approximation that 
we use here is better than the one obtained using a two- 
band low-/c approximation. Indeed, for bilayer graphene 
where both the two-band and the full band polarizations 
were calculated, we see that the \ow-k approximation 
of the two-band model misses the correct high- A: linear 
asymptotics and introduces a large error in the integrals 
which are performed up to the cutoff A. 

Let us now use the linear expression for II and define 
V = gV. Then, since U tot V = 2UV is constant in k, 
the renormalized potential can be written as Vij = gVij, 
where 
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Thus, the large momentum behavior of the renormalized 
potential effectively renormalizes the electron-electron 
coupling g. Let n(g) be an interpolation function rep- 
resenting the phase boundary in the case of no screening 
(the solid line in Fig. [5|. Then, 



n\g) n(g) n 
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is the phase boundary in the screened case. This bound- 
ary is shown by a dashed line in Fig. [5] 

The \ow-k regime of the polarization H(Q a ,k, 0) for 
A£?C-trilayer graphene can be approximated by a con- 
stant w — —l/[67rk F (3], where ft = 400 We will con- 
sider the case where g = 6 and use the critical doping 
k F w 0.0116 (see Fig. [§, which leads to w « -0.0114. 
It is natural to let the transition into the linear regime 
of II occur at the point where nk^ = w, i.e. at fco = 0.1 
for g = 6. The renormalized potential for k < ko now 
becomes 

Vnm(k) V nm (k) - _Jl rgw / k Vnm(k)g(k). 
As a crude approximation, we can let 

Vnm ~ Vnm aVg fcG ^ [#(&)] 

Vnm f k ° f k ° , f 2n g 

vol tt J dp J dp Jo 1 - Airgw/kfap', 0) ' 

where Q = [0,/co] 2 x [0, 2n] is the domain where the 
constant regime of the polarization holds and k = 
y/p 2 + p' 2 — 2pp' cos 0. Thus, the renormalized g for 
g = 6 at the critical fc^ becomes 

where the integration was calculated numerically. At the 
same values of g and k F , the polarization in the linear 
regime yields g « 0.58. Thus, g ~ g, which implies 
that, as a first approximation, we may consider the linear 
approximation of the polarization for all momenta, which 
leads to the phase boundary represented by the black 
dashed line of Fig. [5] 

V. CONCLUSIONS 

In this paper we study the magnetic properties of 
A£?C-trilayer graphene using a tight-binding approach, 
where only the nearest-neighbor hopping parameters are 
taken into account. We include the Coulomb interaction 
and evaluate the exchange energy (Fock term) allowing 
for an unequal filling of the spin-up and spin-down bands. 
Then we calculate numerically the difference in energy 
between paramagnetic and ferromagnetic configurations 
of the system and identify the points of phase transitions 
for fixed values of the interaction parameter g. By re- 
peating the calculations for several values of g, we obtain 
the phase diagram in the electron-electron coupling vs. 
doping plane. As a first step, we did not take into ac- 
count the effects of Coulomb screening. The results are 
shown as the solid line in Fig. [5] 

Although the phase diagram for monolayer, 27 
bilayer, 28 and Al^-trilayer 3 -^ graphene have been 
previously derived, effects of screening have been ne- 
glected until now. Our work represents the first step to 
incorporate these important effects. 



For the unscreened case, at g « 2.1, a comparison 
with unscreened bilayer graphene^ shows that ABC- 
trilayer graphene has a ferromagnetic behavior which is 
approximately 50 times stronger. Furthermore, a simi- 
lar comparison with A£?A-trilayer graphene 30 shows that 
A£?C-trilayer has a ferromagnetic behavior that is ap- 
proximately 300 times stronger. At g « 2.1, monolayer 
graphene shows a paramagnetic behavior at all doping 
levels. In order for phase transitions to be present in 
monolayer graphene, the electron-electron coupling needs 
to exceed g « 5. 27 Fig. [Z] shows that at g = 6, the 
phase transition in A^BC-trilayer graphene is of first or- 
der. This behavior persists for all couplings g < 6. ABA- 
trilayeJ^ and bilayer graphene 28 also exhibits first order 
phase transitions for couplings g < 6. This is in contrast 
to monolayer graphene, where both first order and sec- 
ond order phase transitions take place at given couplings 
gW\ Thus, A£?C-trilayer graphene behaves in a similar 
manner to bilayer and ABA-triiayer graphene, but ex- 
hibits a much stronger ferromagnetic behavior, making 
it easier to experimentally detect ferromagnetism. At 
g = 2.1, the phase transition to ferromagnetism occurs 
at n w 5.5 • 10 -5 . In Si-units the doping level becomes 
n = # s ^A 2 Q 2 /[4tt] = g s g v Q 2 d /[2A] = ng s g v /A, where 
g s = 2 and g v = 2 are the spin and valley degenera- 
cies, respectively, and A w 5.2 • 10 -16 cm 2 is the area of 
the Brillouin zone. Thus, neglecting valley degeneracy, 
ft ~ 2 • 10 11 cm -2 . Note that, by mapping the parameter 
x' of Fig. [I] to x, we see that the critical curve attains a 
minimum at x < 0. Thus, in the ferromagnetic regime at 
g = 6, the energy is always minimized for a configuration 
with two types of charge carriers. This behavior persists 
for all g < 6. 

These conclusions were reached by neglecting Coulomb 
screening. However, due to the diverging density of states 
in A^C-trilayer graphene, screening plays a very impor- 
tant role and must be taken into account. A thorough 
calculation of the polarization bubble in the full-band 
model is beyond the scope of this paper, and will be de- 
ferred to a future publication. 39 Nevertheless, we have 
included screening effects within a simplified model. In 
the case of monolayer and bilayer graphene, the large-Zc 
behavior of the bubble diagrams are linear in fc, with 
the same slope k. Arguing that this linear behavior 
also applies to A£?C-trilayer graphene, and approximat- 
ing the low-fc behavior of the polarization by a constant, 
we found that screening effects can be incorporated via a 
simple renormalization of the electron-electron coupling 
g. Fig.|5]shows that the large momentum behavior of the 
screening leads to a reduced ferromagnetic region in the 
A£?C-trilayer graphene phase diagram. However, ferro- 
magnetism is still approximately 25 times stronger than 
in unscreened bilayer graphene, which means that ABC- 
trilayer remains the material with the strongest ferro- 
magnetic behavior. 

We are aware that next-nearest neighbor hopping pa- 
rameters, like 73 of the SWMc model can be of the same 
order as 71 p°J and that this can have an influence on the 
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low-momentum behavior of the model. This parameter 
has been systematically neglected in studies of ferromag- 
netism in multilayer graphene (see Ref. 28 for bilayer and 
Ref . l30l for AB A-triiayer) . The reason is that, for study- 
ing the effects of other hopping parameters, one needs to 
redefine what is meant by a particle and a hole pocket 
due to the broken rotational symmetry of the dispersion 
around the K-poi nt of the Brillouin zone, resulting from 
the SWMc mo del ESEU Furthermore, this broken symme- 
try leads to more complex integration boundaries, which 
makes the resulting numerical integrations intractable. 

Recently, an intrinsic bandgap of 6 meV was experi- 
mentally observed in suspended A£?C-trilayer graphene, 
and it was argued that it should be driven by 
interactions However, this gap did not appear in most 
of the samples placed on a substrate, which were inves- 
tigated during the same study. Since suspended samples 
are more susceptible to ripples and deformations, it can 
well be that the spatial inversion symmetry was broken 
by strain, resulting in the intrinsic bandgap. Our stud- 
ies should then apply for ABC-tnlayei graphene on a 
substrate, without deformations. Because the dielectric 
constant is larger for samples on a substrate than for 
suspended samples (in vacuum), the coupling constant g 
will be renormalized by a factor e ~ 2.5 for graphene on a 
Si02 wafer. Otherwise, the paramagnetic-ferromagnetic 
phase transition remains unaltered. 

A simplified theoretical model which includes only on- 
site interactions suggests that the difference in bandstruc- 
ture between ABA- and ABC- stacked trilayers should be 
enough to explain the presence of a gap due to antiferro- 
magnetism in ABC samples, while AB A-stacked trilay- 
ers remain ungappedP^I These studies, however, cannot 
explain why the gap arises only in suspended samples. 

Here we include long-range Coulomb interactions and 
investigate also the effect of screening. It is usually 
argued (without further ado) that screening is more 
important in A^BC-trilayer than in the other related 
compunds. Our studies reveal that this is not always 
true, since the polarization is linearly increasing in a 
considerable region, over which one must integrate to 
obtain the exchange energy. This feature is similar in 
monolayer, bilayer, and ABC- trilayer graphene, and 
it is simply a consequence of the linear dispersion at 
intermediate values of /c, which occurs in all the cases. 
Our studies reveal that the low energy approximation 
for the polarization is not always enough to ground 
fast conclusions. Although the final understanding 
about ABC- trilayer graphene has not yet been reached, 
we hope that our work will pave the way to possible 
extensions of the existing models for the investigation of 
ferromagnetism in multi-layer graphene using numerical 
methods. 
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Appendix A: Exchange energy 

The interaction Hamiltonian for 74£?C-trilayer 
graphene is shown in Eq. (J3|. Fourier transforming p n 
and tycr,n leads to 

P M) = J rfrp„(r)e-^ r = £ j dr^r^^e-^ 



7E^,«( k+( i)*^( k )' 



A 



(Al) 



k,cr 



in the discrete limit. Using Eq. (Al) and Fourier trans- 
forming V D , V ND and V 2ND in Eq. going to the 
discrete limit, and subsequently rewriting the resulting 
expression into a matrix form yields 



i ^ Mq)\ 

^7 22 (M-q) P2(-q) fl»(-q)) M p 2 (q) 



(A2) 



where (by omitting the q dependence for brevity) 
M = M t + M r 

V 2Nm 


y2ND q o 



The matrix M t is diagonalized by U t such that U^D t U t 
M t where 



V D 


yND 


\ 


yND 


V D 


yND 





yND 


V D j 

















V2 














with Vl = V D , v 2 = V D - V2V ND and v 3 = V D + 
\/2V ND . Similarly, M r is diagonalized by U r such that 
UjD r U r = M r where 



1 
U r = I -l/y/2 1/V2 | , D r 
1/V2 1/V2J 



'0 

— y2ND 
,0 V 2ND 




~Pl + P3 

-j= [pi/V2- P 2 + p 3 /V2] , 
V2 \ P i/V2 + p2+p 3 /V2j 




Using the above diagonalizations yields 
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-y 



27re 2 1 



2A ^o ^ ^ 



(Pi(-q) /5 2 (-q) p 3 (-q)) 1 - y/2e~« d 



-y 



27re 2 1 



2A^ 



'0 



(P4(-q) P 5 (-q) Pe(-q)) -e 



-2qd 



,0 




-2qd ) y/2 



By defining 



Let 



VM = ^ V 2/3 (q) = — (1=fV2c-^), 



Vk((z) = 0, Vj>/ 6 (</) 



eg 



-( T e- 2 ^), 



the Hamiltonian reduces to the compact form 
1 6 

Hi = ^yyPs(-(l)V s ( q )p s (q). 

Inspection of the operators p s for s = 1, 2, . . . , 6, indicates 
that they are all linear combinations of p n (q) for n = 
1,2,3. Thus, using Eq. (Al) one obtains 



p s (q) = £*t(k + q)x s *(k) 
k 

= £V( k + q)^ f (k + q)x s X(k)$(k), 

k 

where .M(q) is the diagonalizing matrix of the ABC- 
trilayer Hamiltonian, § = (^1,^2,^3), with § n (q) = 
\I/ n (q)/V^4 being a two component dimensionless anni- 
hilation operator working on layer n. The operator & 
contains the band creation operators of the six bands. 
Thus, the six matrices x a are defined as 



Xl/6 



X2/3 



'ti 2 

1 , x 4 

1 2 , 

'l 2 />/2 

TI2 
1 2 /V2/ 



'0 0> 
y/21 2 

,0 0; 



where %5 = Xi- By defining x s = A^k + q^gA^k) the 
A£?C-trilayer interaction Hamiltonian can be written as 

^=^EEE E ^(p-q)x^(p-q, P ) 
x ^(p)y s (q)$t (p' + q)*f „(p' + q,p')<Mp')- 



|N)= J] [^(k)]^-- |0> 

k,zy,cr 



denote a Fock state of the system, where N^ a ^ G {0, 1} 
is the occupancy of electrons in the momentum state k 
of energy band v with spin a. Then, to first order, the 
energy of the system is described by 



E T = <N| : B l : |N), 



where : : denotes normal ordering. Working out the 
expectation value of Hj results in two distinct contri- 
butions. These are the Hartree (direct) and the Fock 
(exchange) contributions. Because of the Jellium back- 
ground the Hartree contribution vanishes and only the 
Fock contribution remains. Thus, 



^ DO 

£ - = ^EE E E^a^P'P') 

p,p' s=l a, (3=1 cr,a 

x X^(P,P')^(P'-P)X^ Q (P',P), 



where 



Aa tatPt a = " (p)*^ ( P )*t ^ (p')d>^ (p') | N ) 

= -^cT,/3,a(pO n cT,a,a(p), 



and tl<j a a(p') &re Fermi occupation functions, which in 
the T —> limit become Heaviside step functions repre- 
senting the pocket configurations shown in Fig. [3j Going 
to the continuum limit reproduces the result shown in Eq. 
(J6|. For further information on the numerical methods 
used to solve the exchange integral, see Ref. [42j 



10 



* richard.olsen.75@gmail.com 

1 K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. 
Zhang, S.V. Dubonos, I.V. Grigorieva, and A. A. Firsov, 
Science 306, 666 (2004). 

2 Keun Soo Kim, Yue Zhao, Houk Jang, Sang Yoon Lee, 
Jong Min Kim, Kwang S. Kim, Jong-Hyun Ahn, Philip 
Kim, Jae- Young Choi, and Byung Hee Hong, Nature 457, 
706 (2009). 

3 Helin Cao, Qingkai Yu, Robert Colby, Deepak Pandey, 
C. S. Park, Jie Lian, Dmitry Zemlyanov, Isaac Childres, 
Vladimir Drachev, Eric A. Stach, Muhammad Hussain, 
Hao Li, Steven S. Pei, and Yong P. Chen, J. Appl. Phys. 
107, 044310 (2010). 

4 Sukang Bae, Hyeongkeun Kim, Youngbin Lee, Xiang- 
fan Xu, Jae-Sung Park, Yi Zheng, Jayakumar Balakrish- 
nan, Tian Lei, Hye Ri Kim, Young II Song, Young-Jin 
Kim, Kwang S. Kim, Barbaras Ozyilmaz, Jong-Hyun Ahn, 
Byung Hee Hong, and Sumio Iijima, Nature Nanotech. 5, 
574 (2010). 

5 J. L. Tedesco, B. L. VanMil, R. L. Myers-Ward, J. M. 
McCrate, S. A. Kitt, P. M. Campbell, G. G. Jernigan, J. 
C. Culbertson, C. R. Eddy, Jr., and D. K. Gaskill, Appl. 
Phys. Lett. 95, 122102 (2009). 

6 Kenjiro K. Gomes, Warren Mar, Wonhee Ko, Francisco 
Guinea, and Hari C. Manoharan, Nature 483, 306 (2012). 

7 Alfonso Reina, Xiaoting Jia, John Ho, Daniel Nezich, 
Hyungbin Son, Vladimir Bulovic, Mildred S. Dresselhaus, 
and Jing Kong, Nano Lett. 9, 30 (2009). 

8 K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, 
M.I. Katsnelson, I.V. Grigorieva, S.V. Dubonos, and A. A. 
Firsov, Nature 438, 197 (2005). 

9 K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal'ko, 
M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, and A. 
K. Geim, Nature Phys. 2, 177 (2006). 

10 M. Koshino and E. McCann, Phys. Rev. B 81, 115315 

(2010) . 

11 Hongki Min and A.H. MacDonald, Phys. Rev. B 77, 
155416 (2008). Hongki Min and A.H. MacDonald, Prog. 
Theor. Phys. Suppl. 176, 227 (2008). 

12 J.C. Slonczewski and P.R. Weiss, Phys. Rev. 109, 272 
(1958). 

13 J.W. McClure, Phys. Rev. 108, 612 (1957). 

14 Thiti Taychatanapat, Kenji Watanabe, Takashi Taniguchi, 
and Pablo Jarillo-Herrero, Nature Phys. 7, 621 (2011). 

15 S. H. Jhang, M. F. Craciun, S. Schmidmeier, S. Toku- 
mitsu, S. Russo, M. Yamamoto, Y. Skourski, J. Wosnitza, 
S. Tarucha, J. Eroms, and C. Strunk, Phys. Rev. B 84, 
161408(R) (2011). 

16 R. Ma, L. Sheng, M. Liu, and D.N. Sheng, Phys. Rev. B 
86, 115414 (2012). 

17 S. Bala Kumar and Jing Guo, Appl. Phys. Lett. 100, 
163102 (2012). 

18 Chun Hung Lui, Zhiqiang Li, Kin Fai Mak, Emmanuele 
Cappelluti, and Tony F. Heinz, Nature Physics, 7, 944 

(2011) . 

19 M. Koshino and E. McCann, Phys. Rev. B 80, 165409 
(2009). 

20 Jia-An Yan, W. Y. Ruan, and M. Y. Chou, Phys. Rev. B 
77, 125401 (2008). 



21 Zhiqiang Li, Chun Hung Lui, Emmanuele Cappelluti, Lara 
Benfatto, Kin Fai Mak, G. Larry Carr, Jie Shan, and Tony 
F. Heinz, Phys. Rev. Lett. 108, 156801 (2012). 

22 A. Kumar, W. Escoffier, J.M. Poumirol, C. Faugeras, D.P. 
Arovas, M.M. Fogler, F. Guinea, S. Roche, M. Goiran, and 
B. Raquet, Phys. Rev. Lett. 107, 126806 (2011). 

23 Fan Zhang, Dagim Tilahun, and A. H. MacDonald, Phys. 
Rev. B 85, 165139 (2012). 

24 Kin Fai Mak, Jie Shan, and Tony F. Heinz, Phys. Rev. 
Lett. 104, 176404 (2010). 

25 Jia-An Yan, W. Y. Ruan, and M. Y. Chou, Phys. Rev. B 
83, 245418 (2011). 

26 J. Borysiuk, J. Soltys, and J. Piechota, Nature 109, 093523 
(2011). 

27 N.M.R. Peres, F. Guinea, and A.H. Castro Neto, Phys. 
Rev. B 72, 174406 (2005). 

28 J. Nilsson, A.H. Castro Neto, N.M.R. Peres, and F. 
Guinea, Phys. Rev. B 73, 214418 (2006). 

29 J. Martin, N. Akerman, G. Ulbricht, T. Lohmann, J. H. 
Smet, K. von Klitzing, and A. Yacoby, Nature Physics 4, 
144 (2008). 

30 Ralph van Gelderen, Lih-King Lim, and C. Morais Smith, 
Phys. Rev. B 84, 155446 (2011). 

31 Mikito Koshino, Phys. Rev. B 81, 125304 (2010). 

32 We are using the letter a to denote both, the lattice con- 
stant and the valley degree of freedom, as it is usually done 
in the literature. 

33 Masatake Mori and Masaaki Sugihara, J. Comput. Appl. 
Math. 127, 287-296 (2001). 

34 Michael G. Duffy, SIAM J. Numer. Anal. 19, 1260-1262 
(1982). 

35 William H. Press, Saul A. Teukolsky, William T. Vet- 
terling, and Brian P. Flannery, Numerical Recipes, Cam- 
bridge, 2007. 

36 Henrik Bruus and Karsten Flensberg, Many-Body Quan- 
tum Theory in Condensed Matter Physics, Oxford, 2004. 

37 O. V. Gamayun, Phys. Rev. B 84, 085112 (2011). 

38 There is a factor of l/(27r) difference from the article of 
Gamayun, due to differing conventions of the Fourier trans- 
form. 

39 Ralph van Gelderen, Richard Olsen, and C. Morais Smith, 
unpublished (2012). 

40 F. Zhang, B. Sahu, H. Min, and A. H. MacDonald, Phys. 
Rev. B 82, 035409 (2010). 

41 Ralph van Gelderen and C. Morais Smith, Phys. Rev. B 
81, 125435 (2010). 

42 Richard Olsen, Ferromagnetism in ABC -trilayer graphene. 
Mast er's thesis (2012), 



http://web.science.uu.nl/ITF/Teaching/2012/Richard 



%20Olsen.pdf (includes updates and corrections). 
W. Bao, L. Jing, J. Velasco Jr, Y. Lee, G. Liu, D.Tran, B. 
Standley, M. Aykol, S.B. Cronin, D. Smirnov, M. Koshino, 
E. McCann, M. Bockrath, and C.N. Lau, Nat. Phys. 7, 948 
(2011). 

D.-H Xu, J. Yuan, Z.-J Yao, Y. Zhou, J.-H Gao, and F.-C 
Zhang, Arxiv: 1207.5287 (2012). 



